!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate thrust produced by the nozzle.
! ref. to chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-07-07
! revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : Thrust
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   pp4   pressure on the wall
!   yy4   radical location of the point
!   Mtemp Mach number of the point
!   pw    front wall point oressure on the wall
!   yw    front wall point radical location on the wall
!   Fm    mass flow momentum thrust
!   Fp    wall pressure thrust
!   dF    thrust act through the point
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine Thrust(pp4,yy4,MM4,pw,yw)
   use VariableDef
   implicit none
!   Given variables
   real*8::pp4,yy4,MM4,pw,yw
!   temp variables
   real*8::Mtemp,eetemp,Fm,Fp,Mp,F1,F2,dF
   real*8::temp1,temp2,temp3,temp4
!   calculate variavbles
!   real*8::F1d,F,etaF,etaI,Isp,Isp1d
   
   Mtemp=MM4
   eetemp=(yy4/yt)**2
   temp1=(g-1)/2
   temp2=(g+1)/2
   temp3=temp2/(g-1)
   temp4=(1/temp2)**temp3

!   one dimensional thrust
   do
      Fm=Mtemp*eetemp-temp4*(1+temp1*Mtemp**2)**temp3
      Fp=eetemp-temp4*temp2*Mtemp*(1+temp1*Mtemp**2)**(temp3-1)
      Mp=Mtemp-Fm/Fp
      if (abs(Mp-Mtemp) < 1e-6) exit
      Mtemp=Mp
   end do
   F1d=Fs*(1+g*Mp**2)/sqrt(2*(g+1)*Mp**2*(1+temp1*Mp**2))-p0*pi*yy4**2   ! 1-d thrust
   
!   two dimensional thrust
   F1=pi*(pw-p0)*yw
   F2=pi*(pp4-p0)*yy4
   dF=(F1+F2)*(yy4-yw)
   F=F+dF   !thrust add
   
   etaF=F/F1d
   etaI=etaF/CDm
   Isp=F/mdot   ! specific impulse
   Isp1d=F1d/mdot1d   ! 1-d specific impulse
end subroutine Thrust
